include("operators.jl")

function MeshReaderGambit2D(filename)
    """
    Read in basic grid information to build grid
    gambit(Fluent, Inc) *.neu format is assumed
    """
    # Open the .neu format meshfile
    file = open(filename)

    # Neglect the head of .neu file
    for i = 1:6
        line = readline(file)
    end

    # Get number of nodes, cells
    dim = parse.(Int, split(strip(readline(file)),r"\s+"))
    Nv, K = dim

    # Skip two lines
    for i = 1:2
        line = readline(file)
    end

    # Read coordinates of the nodes
    Vx = [ ]
    Vy = [ ]
    coordinates = zeros(Float64,Nv,2)
    for i in 1:Nv
        line = parse.(Float64,split(strip(readline(file), ' '),r"\s+"))
        coordinates[i,1] = line[2]
        coordinates[i,2] = line[3]
    end

    # Skip two lines
    for i = 1:2
        line = readline(file)
    end

    # Read elements to nodes connectivity
    EtoV = zeros(Int, K, 3)
    for i = 1:K
        line = parse.(Int,split(strip(readline(file), ' '),r"\s+"))
        EtoV[i,1] = line[4]; EtoV[i,2] = line[5]; EtoV[i,3] = line[6]
    end

    # Skip the material property section
    for i = 1:4
        line = readline(file)
    end
    test = intersect(line, ["ENDOFSECTION"])
    while isempty(test)
        line = split(strip(readline(file), ' '),r"\s+")
        test = intersect(line, ["ENDOFSECTION"])
    end

    line = readline(file); line = readline(file)

    # Define boundary type
    BCType = zeros(Int, K,3)
    In = 1; Out = 2; Wall = 3; Far = 4; Cyl = 5; Dirichlet = 6; Neuman = 7; Slip = 8;Up = 9; Low=10;

    while ~isempty(line)
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["In"])) bcflag = In end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Out"])) bcflag = Out end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Wall"])) bcflag = Wall end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Far"])) bcflag = Far end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Cyl"])) bcflag = Cyl end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Dirichlet"])) bcflag = Dirichlet end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Neuman"])) bcflag = Neuman end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Slip"])) bcflag = Slip end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Up"])) bcflag = Up end
        if ~isempty(intersect(split(strip(line, ' '),r"\s+"), ["Low"])) bcflag = Low end

        line = split(strip(readline(file), ' '),r"\s+")
        while isempty(intersect(line, ["ENDOFSECTION"]))
            boundary_id = parse.(Int, line)
            BCType[boundary_id[1], boundary_id[3]] = bcflag
            line = split(strip(readline(file), ' '),r"\s+")
        end
        line = readline(file); line = readline(file)
    end

    close(file)
    Nv, K, coordinates, EtoV, BCType
end

function Connectivity2D(EToV,Nfaces,K,Nv)
    """
    Purpose: triangle face connect algorithm due to Toby Isaac
    """
    # Create list of all faces 1, then 2, & 3
    fnodes = [EToV[:,[1,2]]; EToV[:,[2,3]];EToV[:,[3,1]]]
    fnodes = sort(fnodes, dims = 2) .- 1

    # Set up default element to element and Element to faces connectivity
    EToE= collect(1:K)*ones(Int,1,Nfaces)
    EToF= ones(Int,K,1)*(1:Nfaces)'

    # Uniquely number each set of three faces by their node numbers
    id = fnodes[:,1]*Nv + fnodes[:,2] .+ 1
    spNodeToNode = [id collect(1:Nfaces*K) EToE[:] EToF[:]]
    sorted = sortslices(spNodeToNode, dims=1)

    # Find matches in the sorted face list
    indices = findall(sorted[1:end-1,1].== sorted[2:end,1])

    # make links reflexive
    matchL = [sorted[indices,:]  ; sorted[indices .+ 1,:]]
    matchR = [sorted[indices .+ 1,:]; sorted[indices,:]]

    # insert matches
    EToE[matchL[:,2]] = matchR[:,3]
    EToF[matchL[:,2]] = matchR[:,4]

    EToE, EToF
end

function Lift2D(N,Np,Nfaces,Nfp,Fmask,r,s,V)
    """
    Purpose  : Compute surface to volume lift term for DG formulation
    """
    Emat = zeros(Np, Nfaces*Nfp)

    # face 1
    faceR = r[Fmask[:,1]]
    V1D = Vandermonde1D(N, faceR)
    massEdge1 = inv(V1D*V1D')
    Emat[Fmask[:,1],1:Nfp] = massEdge1

    # face 2
    faceR = r[Fmask[:,2]]
    V1D = Vandermonde1D(N, faceR)
    massEdge2 = inv(V1D*V1D')
    Emat[Fmask[:,2],Nfp+1:2*Nfp] = massEdge2

    # face 3
    faceS = s[Fmask[:,3]]
    V1D = Vandermonde1D(N, faceS)
    massEdge3 = inv(V1D*V1D')
    Emat[Fmask[:,3],2*Nfp+1:3*Nfp] = massEdge3

    # inv(mass matrix)*\I_n (L_i,L_j)_{edge_n}
    LIFT = V*(V'*Emat)

    LIFT
end

mutable struct Geometry
    """
    A struct to handle geomatric informations.
    """

    # mesh informations
    K::Int64                  # number of elements
    Nv::Int64                 # number of vortices in grid
    Nfp::Int64                # number of Flux Points at each face
    NFaces::Int64             # number of faces in each element

    Vs::Array{Float64,2}      # locations of Element Vortices


    # connectivity informations, all of integer type
    EtoV::Array{Int64,2}      # array of Element to Vortices
    EtoF::Array{Int64,2}      # array of Element to Faces
    EtoE::Array{Int64,2}      # array of Element to Elements
    BCType::Array{Int64,2}    # array of boundary condition types

    vmapM::Array{Int64,1}     # vector of global nodal numbers at face for interior values
    vmapP::Array{Int64,1}     # vecotr of global nodal numbers at face for exterior values
    vmapB::Array{Int64,1}     # vecotr of global nodal numbers at face for boundary values

    vmapI::Array{Int64,1}
    vmapO::Array{Int64,1}
    vmapU::Array{Int64,1}
    vmapL::Array{Int64,1}

    mapB::Array{Int64,1}
    mapP::Array{Int64,1}

    mapI::Array{Int64,1}
    mapO::Array{Int64,1}
    mapU::Array{Int64,1}
    mapL::Array{Int64,1}

    # get solution points and flux points
    SPs::Array{Float64,3}     # locations of Solution Points(SPs)
    FPs::Array{Float64,3}     # locations of Flux Points at(FPs)
    FPidx::Array{Int64,2}     # Flux Points' coordinates

    # Operators and local matrics
    rx::Array{Float64,2}
    sx::Array{Float64,2}
    ry::Array{Float64,2}
    sy::Array{Float64,2}
    Jv::Array{Float64,2}     # Volume Jacobian, J
    Js::Array{Float64,2}     # Jacobian at surface nodes

    Normal::Array{Float64,3} # Normal vector at FPs
    Fscale::Array{Float64,2} # Fscale(i,k) is the ratio of surface to volume Jacobian
    LIFT::Array{Float64,2}   # create surface integration term

    function Geometry(mesh,op)
        self = new()

        # get basic mesh informations
        Nfp = op.P + 1; self.Nfp = Nfp
        self.NFaces = 3; Nfaces = self.NFaces
        self.Nv,self.K,self.Vs,self.EtoV,self.BCType = MeshReaderGambit2D(mesh)
        self.EtoE,self.EtoF = Connectivity2D(self.EtoV,self.NFaces,self.K,self.Nv)


        # get solution points(SPs) distribution
        Vx = self.Vs[:,1]; Vy = self.Vs[:,2]
        va = self.EtoV[:,1]'; vb = self.EtoV[:,2]'; vc = self.EtoV[:,3]'

        self.SPs = zeros(Float64,op.Np,self.K,2)
        self.SPs[:,:,1] = 0.5*(-(op.r+op.s)*Vx[va] + (1 .+ op.r)*Vx[vb] + (1 .+ op.s)*Vx[vc])
        self.SPs[:,:,2] = 0.5*(-(op.r+op.s)*Vy[va] + (1 .+ op.r)*Vy[vb] + (1 .+ op.s)*Vy[vc])


        # get Flux Points(FPs) distribution
        NODETOL = 1e-12
        Face₁ = findall(abs.(op.s[:] .+ 1) .< NODETOL)
        Face₂ = findall(abs.(op.r[:] .+ op.s[:]) .< NODETOL)
        Face₃ = findall(abs.(op.r[:] .+ 1) .< NODETOL)

        SPx = self.SPs[:,:,1]; SPy = self.SPs[:,:,2]
        self.FPidx = [Face₁ Face₂ Face₃]

        self.FPs = zeros(Float64,self.Nfp*3,self.K,2)
        self.FPs[:,:,1] = SPx[self.FPidx[:],:]
        self.FPs[:,:,2] = SPy[self.FPidx[:],:]


        # calculate geomatric metrics
        xr = op.Dr*SPx; xs = op.Ds*SPx
        yr = op.Dr*SPy; ys = op.Ds*SPy

        self.Jv = -xs.*yr + xr.*ys
        self.rx =  ys./self.Jv; self.sx = -yr./self.Jv
        self.ry = -xs./self.Jv; self.sy =  xr./self.Jv


        # initial the normal vector array
        nx = zeros(Nfp*3,self.K); ny = zeros(Nfp*3,self.K)

        # interpolate geometric factors to face nodes
        Fxr = xr[self.FPidx[:], :]; Fxs = xs[self.FPidx[:], :]
        Fyr = yr[self.FPidx[:], :]; Fys = ys[self.FPidx[:], :]

        # set Ids
        Id₁ = 1:Nfp; Id₂ = Nfp .+ 1:2*Nfp; Id₃ = 2*Nfp .+ 1:3*Nfp

        # face 1
        nx[Id₁, :] =  Fyr[Id₁, :]
        ny[Id₁, :] = -Fxr[Id₁, :]

        # face 2
        nx[Id₂, :] =  Fys[Id₂, :] - Fyr[Id₂, :]
        ny[Id₂, :] = -Fxs[Id₂, :] + Fxr[Id₂, :]

        # face 3
        nx[Id₃, :] = -Fys[Id₃, :]
        ny[Id₃, :] =  Fxs[Id₃, :]

        # normalize
        self.Js = sqrt.(nx.*nx + ny.*ny)
        nx = nx./self.Js; ny = ny./self.Js

        # get the normal vectors
        self.Normal = zeros(Nfp*3,self.K,2)
        self.Normal[:,:,1] = nx; self.Normal[:,:,2] = ny

        # get face fscaling
        self.Fscale = self.Js ./ self.Jv[self.FPidx[:],:]

        # get surface integration term
        self.LIFT = Lift2D(op.P,op.Np,Nfaces,Nfp,self.FPidx,op.r,op.s,op.V)


        # build up mesh, get left and right side status
        Np = op.Np; K = self.K

        ids   = reshape(1:K*Np, Np, K)
        vmapM = zeros(Int, Nfp, Nfaces, K)
        vmapP = zeros(Int, Nfp, Nfaces, K)
        mapM  = collect(Int, 1:K*Nfp*Nfaces)
        mapP  = reshape(collect(Int, 1:K*Nfp*Nfaces),Nfp, Nfaces, K)

        # find index of face nodes with respect to volume node ordering
        for k=1:K
            for f=1:Nfaces
                vmapM[:,f,k] = ids[self.FPidx[:,f], k]
            end
        end

        one = ones(1, Nfp)
        for k₁=1:K
          for f₁=1:Nfaces
            # find neighbor
            k2 = self.EtoE[k₁,f₁];  f₂ = self.EtoF[k₁,f₁]

            # reference length of edge
            v1 = self.EtoV[k₁,f₁]; v2 = self.EtoV[k₁, 1 + mod(f₁,Nfaces)]
            ref = sqrt.((Vx[v1]-Vx[v2]).^2 + (Vy[v1]-Vy[v2]).^2)

            # find volume node numbers of left and right nodes
            vidM = vmapM[:,f₁,k₁]; vidP = vmapM[:, f₂,k2]
            x1 = SPx[vidM]; y1 = SPy[vidM]; x2 = SPx[vidP]; y2 = SPy[vidP]
            x1 = x1*one;  y1 = y1*one
            x2 = x2*one;  y2 = y2*one

            # compute distance matrix
            D = (x1 -x2').^2 + (y1-y2').^2

            index = findall(sqrt.(abs.(D)) .< NODETOL*ref)
            idM = [index[i][1] for i = 1:length(index)]
            idP = [index[i][2] for i = 1:length(index)]

            vmapP[idM,f₁,k₁] = vidP[idP]
            mapP[idM,f₁,k₁] = idP .+ ((f₂-1)*Nfp+(k2-1)*Nfaces*Nfp)
          end
        end

        # reshape vmapM and vmapP to be vectors and create boundary node list
        self.vmapP = vmapP[:]
        self.vmapM = vmapM[:]

        self.mapP = mapP[:]
        self.mapB = findall(vmapP[:] .== vmapM[:])
        self.vmapB = vmapM[self.mapB]


        #Create label of face nodes with boundary types from BCType
        bct    = self.BCType'
        bnodes = ones(Nfp, 1)*bct[:]'
        bnodes = bnodes[:]

        #find location of boundary nodes in face and volume node lists
        self.mapI = findall(bnodes .== 1); self.vmapI = vmapM[self.mapI]
        self.mapO = findall(bnodes .== 2); self.vmapO = vmapM[self.mapO]

        self.mapU = findall(bnodes .== 9);  self.vmapU = vmapM[self.mapU]
        self.mapL = findall(bnodes .== 10); self.vmapL = vmapM[self.mapL]
        return self
    end
end
